# -------------------------------------------------------------------------
# Replication
# Author: Alison Pei
# Date created; May. 9th 2022
# -------------------------------------------------------------------------

#Packages to install -----------------------------------------------------------
# install.packages("did")

set.seed(123)
# Load packages -----------------------------------------------------------
library(tidyverse)
library(AER)
library(haven)
library(lubridate)
library(dplyr)
library(readr) 
library(MASS)
library(rticles)


# Import and Clean data --------------------------------------------------------------
setwd("C:/Users/msj22/Dropbox/research/wrongful_discharge/replication_May_2022/intermediate_dtas/")

data <- read_dta("osha_nsc2_for_cs_estimator.dta")
sample <- subset(data, (state_run_office!=1 & year>=1979) | (state_run_office==1 & year>=1992) )

# Regressions --------------------------------------------------------------
# dependent variable is workplace injury 
#1 reghdfe ln_acc_rateY_matt  wdlapY $weight if sector == 1 & $sample, absorb($fes) vce(cluster state_code) 
library(did)

out1 <- att_gt(yname="ln_acc_rateY_matt",
               tname = "year",
               idname = "state_code",
               gname = "first_year_full_wdlapY",
               weightsname = "sh_emp1979",
               clustervars = "state_code",
               xformla = ~1,
               data=sample)

#summary(out1)

#  Simple Aggregation
#  Return a weighted average of all group-time average treatment effects
#  with weights proportional to the group size. 
agg.simple1 <- aggte(out1, type = "simple")
summary(agg.simple1)


agg.es <- aggte(out1, type = "dynamic")
ggdid(agg.es)



#2reghdfe ln_acc_rateY_matt  wdlapY $weight if sector == 1 & $sample, absorb($fes year_div ) vce(cluster state_code) 

out2 <- att_gt(yname="ln_acc_rateY_matt",
               tname = "year",
               idname = "state_code",
               gname = "first_year_wdlapY",
               weightsname = "sh_emp1979",
               clustervars = "state_code",
               xformla = ~year_div,
               data=sample)

summary(out2)

agg.simple2 <- aggte(out2, type = "simple", na.rm = TRUE)
summary(agg.simple2)

agg.es2 <- aggte(out2, type = "dynamic", na.rm=TRUE)
ggdid(agg.es2)


#3 	reghdfe ln_acc_rateY_matt  wdlapY $weight if sector == 1 & $sample, absorb($fes year_div c.year#i.state_code) vce(cluster state_code) 
out3 <- att_gt(yname="ln_acc_rateY_matt",
               tname = "year",
               idname = "state_code",
               gname = "first_year_wdlapY",
               weightsname = "sh_emp1979",
               clustervars = "state_code",
               xformla = ~ year:state_code,
               data=sample)

#summary(out3)

agg.simple3 <- aggte(out3, type = "simple", na.rm = TRUE)
summary(agg.simple3)

agg.es3 <- aggte(out3, type = "dynamic", na.rm=TRUE)
ggdid(agg.es3)


